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Vertex-Transport  Element  Simulation  of  a 
Confined  Mixing  Layer 

by 
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Kent,  Washington  98032 

ABSTRAC1 

An  improved  vortex-transport  element  numerical  scheme  has  been  employed  to 
simulate  the  spatial  evolution  of  a  two-dimensional  planar,  constant  density 
mixing  layer  formed  behind  a  splitter  plate  configuration.  This  scheme  uses 
vortex  elements  to  discretize  the  region  of  high  vorticity,  and  transport 
elements  to  discretize  the  region  of  high  scalar  gradients.  This  improved 
vortex  method  offers  a  method  of  accurately  resolving  the  areas  of  high  strain 
by  redistributing  the  neighboring  elements  within  those  regions. 

The  resulting  scheme  has  been  used  to  simulate  the  convective  transport  of 
a  conserved  scalar  variable  in  the  shear  layer,  and  the  results  show  favorable 
agreement  with  experimental  data  for  the  mean  and  the  rms  values  of  the  scalar 
quantity.  The  effects  of  harmonic  forcing  are  shown  to  accelerate  the 
mechanism  of  the  rollup  and  pairing  within  the  domain,  which  results  in  an 
increase  in  the  convolution  of  the  surfaces.  This  would  be  of  crucial  interest 
to  the  problem  of  active  combustion  control  in  a  reacting  mixing  layer 
simulation. 


1 


1.  INTRODUCTION 


Nonpremixed  turbulent  reacting  flows  have  been  the  subject  of  extensive 
experimental  and  theoretical  investigations  during  recent  years  [for  reviews, 
see  Bilger  (1976,  1980)  and  Givi  (1984)).  In  most  of  the  theoretical  work, 
turbulence  models  are  used  to  close  a  system  of  averaged  transport  equations, 
which  describes  the  statistical  behavior  of  the  aerothermodynam  ical  variables. 
Moment  methods  (Donaldson  and  Varma,  1976),  eddy  break-up  and  mixing-controlled 
models  (Givi  et  al.,  1984a),  flame  sheet  approximation  (Bilger,  1980),  assumed 
probability  density  function  (PDF)  shape  methods  (Lockwood  and  Naguib,  1975), 
solutions  based  on  modeled  joint  PDF  of  scalar  quantities  (Givi  et  al.,  1984b', 
Nguyen  and  Pope,  1984)  and  based  on  modeled  joint  PDF  of  scalar  and  velocity 
(Pope  and  Correa,  1987)  are  examples  in  which  turbulence  modeling  have  been 
used  for  the  closure  of  equations  governing  the  statistical  quantities.  Much 
effort  has  gone  into  constructing  accurate  models  and  in  obtaining  results  that 
are  in  agreement  with  experimental  measurements.  However,  in  complex  systems, 
modeling  is  difficult  because  of  our  lack  of  knowledge  on  the  detailed  dynamics 
of  the  flow.  Furthermore,  since  most  of  the  interesting  dynamical  behavior  of 
the  flow  is  modeled  a  priori,  such  features  are  not  exhibited  from  the  results 
of  numerical  computations  based  on  turbulence  models,  and  thus  cannot  advance 
our  understanding  of  turbulent  combustion. 

The  progress  in  numerical  methods  and  the  availability  of  supercomputers 
have  had  a  major  impact  on  turbulence  research.  Improved  accuracy  of  the 
numerics  and  increased  storage  and  computational  speed  have  made  it  possible  to 
solve  the  appropriate  transport  equations  governing  turbulent  combustion 
directly  without  the  need  for  modeling  over  some  limited  parameter  range.  Such 
nearly  "model-free"  simulations,  in  comparison  with  calculations  utilizing 
turbulence  models,  have  the  advantage  that  the  dominant  physics  of  the  problem 
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is  not  modeled  a  priori  but  is  recovered  directly  from  the  computed  results. 
These  results  can  be  used  to  understand  many  important  mechanisms  of  turbulent 
transport  and  its  direct  influence  on  chemical  reactions.  Furthermore,  since 
the  instantaneous  behavior  of  the  variables  are  known  at  all  points  and  at  all 
times,  accurate  simulations  offer  a  good  method  of  probing  the  flow  when  exper¬ 
imental  techniques  may  fail.  There  are,  however,  some  limitations  on  the  range 
of  turbulent  scales  that  can  be  resolved  accurately  by  model-free  simulations. 
Therefore,  there  is  a  need  to  validate  the  results  of  the  simulations  by  a 
direct  comparison  with  experimental  measurements.  With  such  validations,  ab 
initio  predictions  can  ultimately  be  a  reality. 

Numerical  methods  have  been  used  in  a  variety  of  forms  for  the  simulation 
of  turbulent  flows  in  complex  configurations.  A  recent  survey  can  be  found  in 
review  articles  (Oran  and  Boris,  1987).  In  reacting  flow,  three  approaches  are 
used:  (1)  finite  difference  methods,  (2)  spectral  methods,  and  (3)  vortex 

methods.  In  the  first  approach,  the  variables  are  defined  on  a  grid  and  the 
transport  equations  are  approximated  by  discretizing  the  derivatives  on  the 
grid  nodes.  Examples  of  this  approach  can  be  found  in  the  work  of  Corcos  and 
Sherman  (1984),  who  used  a  projection  method  to  study  the  temporal  evolution  of 
a  periodic  shear  layer,  and  in  Grinstein  et  al.  (1985),  who  used  a  flux- 
corrected  transport  scheme  to  simulate  the  development  of  coherent  structures 
in  a  two-dimensional  spatially  evolving  shear  layer  and  examined  their  effect 
on  mixing. 

In  spectral  methods,  the  variables  are  expanded  in  a  series  of  harmonic 
functions  that  satisfy  the  differential  equations  on  a  number  of  collocation 
points.  Riley  et  al.  (1986)  used  a  pse  udo- spe  ctral  scheme  to  study  a  three- 
dimensional  tem por al ly- ev olv ing  reacting  mixing  layer  assuming  a  constant 
reaction  rate,  constant  density,  and  no  heat  release.  McMurtry  et  al.  (1986) 
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considered  the  effects  of  the  chemical  heat  release  on  the  flnid  dynamics  of  a 
two-dimensional  mixing  layer  for  a  constant  reaction  rate.  The  interplay 
between  flnid  dynamics  and  the  chemical  reaction  is  investigated  under  these 
conditions.  Givi  et  al.  (1987)  used  the  same  method  to  compute  a  two- 
dimensional  mixing  layer  with  an  Arrhenius  chemical  reaction  and  constant 
density  to  assess  the  effects  of  large  coherent  structure  on  the  local 
extinction  of  the  flame.  Extension  to  spatially-growing  layers  was  initiated 
by  Givi  and  Jou  (1988)  using  a  hybrid  p  seudo- spe  c  tr  al  second  order  finite 
difference  scheme.  In  all  cases,  the  Reynolds  number  was  kept  at  small  values, 
0(100),  limited  by  the  grid  resolution  and  the  number  of  harmonic  modes. 

In  the  third  approach,  vortex  methods  are  used.  These  schemes  are  grid 
free,  the  transport  of  the  variables  takes  place  in  a  Lagrangian  form,  and  the 
solution  is  not  restricted  by  the  geometry  of  the  confinement.  Therefore,  they 
can  provide  accurate  simulations  for  high  Reynolds  number,  spatially  growing 
flows.  Moreover,  vortex  methods  optimize  the  computational  efforts  by  distri¬ 
buting  computational  elements  around  regions  of  high  vorticity.  The  applica¬ 
tion  of  the  method  in  thin  premixed  flame  calculations  with  a  finite  density 
jump  has  been  reported  by  Ghoniem  et  al.  (1982)  and  Sethian  (1984),  among 
others.  In  these  calculations,  the  vortex  method  was  employed  to  compute  the 
flow  field,  and  the  dynamic  effect  of  combustion  was  represented  by  the 
propagation  of  a  thin  interface  at  the  laminar  burning  vel ocity  acting  as  a 
volumetric  source. 

Vortex  methods  were  also  used  in  simulating  diffusion  flames  in  connection 
with  a  finite-difference  approach  for  the  treatment  of  the  scalar  variables. 
Ashurst  and  Barr  (1982)  used  the  vortex  method  to  compute  the  hydrodynamic 
field  and  an  Eulerian  flux-corrected  transport  algorithm  to  compute  the 
diffusion  and  convection  of  a  conserved  Sbvab-Zeldov  ich  scalar  approximating 
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the  shape  and  convolution  of  the  flame  in  the  limit  of  infinitely  fast  chemical 
reaction.  Lin  and  Pratt  (1987)  used  the  random  vortex  method  to  simulate  the 
large-scale  motion  and  a  Monte-Carlo  method  to  calculate  the  time-dependent 
probability  density  function  of  the  scalar  quantities  for  both  gaseous  and 
aqueous  mixing  layers.  The  PDF  transport  equation,  however,  required  a  closure 
model  for  the  molecular  mixing  term.  The  application  of  Lagrangian  methods  for 
reactive  flow  was  initiated  by  Ghoniem  and  Givi  (1987),  who  used  vortex-scalar 
methods  for  the  simulations  of  a  reacting  mixing  layer  at  low  heat  release.  In 
that  work,  vortex  methods  were  used  for  hydrodynamic  simulations,  whereas 
scalar  methods  were  used  to  simulate  the  scalar  (temperature  and  concentration) 
variables.  The  results  of  that  work,  although  very  encouraging,  indicated  that 
one  would  need  an  excessive  number  of  scalar  elements  to  simulate  the  scalar 
variables  accurately.  Further,  the  results  indicated  that  since  the  elements 
(both  the  vortex  and  the  scalar)  tend  to  evacuate  regions  of  high  local  strain, 
better  numerical  methods  are  required  to  simulate  such  flows  more  accurately. 

From  this  short  review,  it  is  clear  that  numerical  simulations  have  played 
an  important  role  in  elucidating  the  physics  of  turbulent  reacting  flows,  and 
that  there  is  a  continuing  need  for  more  direct  simulations  in  order  to  better 
explain  some  of  the  interesting  physical  phenomena  that  have  been  observed  in 
laboratory  experiments. 

In  this  work,  we  extend  the  vortex  methods  to  study  a  nonreacting 
spatially  developing  mixing  layer.  This  is  an  essential  step  before  dealing 
with  the  problem  of  chemical  reactions.  The  vortex  transport  element  method  is 
developed  to  treat  both  the  hydrodynamic  and  the  scalar  field  in  a  Lagrangian 
sense.  The  fact  that  a  chemical  reaction  is  truly  a  Lagrangian  process,  i.e., 
it  occurs  when  the  particles  (or  macroscopic  elements)  interact  as  they  flow, 
motivate  the  implementation  of  Lagrangian  methods  for  simulations  of  high 
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Reynolds  number  reacting  flows.  The  method  is  capable  of  handling  a  wide 
variety  of  initial  and  boundary  conditions  and  is  not  limited  to  simple  flow 
boundaries.  In  this  report,  we  concentrate  on  the  formulation  of  the  model  and 
the  numerical  schemes,  and  present  some  preliminary  validation  studies  and 
interpretations  of  the  results. 

In  the  next  section,  the  geometrical  configuration  of  a  spatially-evolving 
mixing  layer  is  presented.  In  this  section,  the  formulation  of  the  vortex 
method  is  given,  and  some  sample  results  of  the  vorticity  field  are  presented. 
These  results  are  used  to  express  the  shortcomings  associated  with  vortex 
methods.  In  Section  3,  the  improved  vortex  methods  are  explained  together  with 
a  description  of  the  transport-element  methods.  This  report  ends  in  Section  4, 
where  summary  and  conclusions  are  given. 

2 .  VORTEX  METHODS 

Consider  a  two-dimensional  flow  behind  a  splitter-plate  configuration,  as 
shown  in  Figure  1.  In  this  figure,  the  velocity  on  the  top  stream  is  equal  to 
U1  and  the  velocity  of  the  low-speed  bottom  stream  is  U2.  The  flow  is  assumed 
incompressible,  and  the  effects  of  boundary  layer  near  the  top  and  bottom  walls 
are  neglected. 

Vortex  methods  are  usually  very  appropriate  to  be  used  for  the  simulations 
of  this  type  of  flow  since,  at  least  initially,  the  vorticies  are  concentrated 
within  a  small  region  of  the  flow.  The  formulation  of  this  scheme,  which  has 
been  explained  in  detail  in  the  literature,  will  be  only  briefly  summarized 
here  for  the  purpose  of  continuity  with  the  next  section. 

The  unsteady  hydrodynamic  transport  equations  are  nondimensional ized  with 
respect  to  the  channel  height  H  and  the  velocity  difference  across  the  layer 
Delta  U  =  U1-U2,  Under  these  conditions,  the  normalized  equations  are 
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(1) 


3«d  _  1  2 

T7  +  n  •  Vu>  = 


at  - 


5  5  =  -«o 


Re 


V  <D 


(2) 


where  w  is  the  vorticity,  5  is  the  stream  function  and  U  and  V  are  the  x  and  y 
components  of  the  vorticity,  i.e.. 


u  =  f£ 

3y 


V=  -ii 

3x 


(3) 


(4) 


To  implement  the  vortex  method,  the  vorticity  field  is  represented  by  a 
finite  number  of  vortex  elements  of  finite  cores: 


w(jc.t)  =  J  Ti/62  f(x  -  xi) 


(5) 


where  Ti  is  the  circulation  of  a  vortex  element  and  6  is  the  core  radius,  while 
xi  is  the  center  of  the  element.  f  represents  the  vorticity  distribution 
associated  with  a  vortex  element,  or  the  core  function  (Chorin,  1973',  Bald, 
1979*,  and  Beale  and  Majda,  1982).  The  velocity  field  is  obtained  by  solving 
Equation  (2)  using  the  discrete  vorticity  distribution: 

fi  “  J  Ti  K(jt  -  ji)  +  Up  (6) 

where  K(;r)  =  -(y,-x)/ is  the  kernel  of  the  Poisson  equation,  K(x)  =  J  r 

dr  is  the  circulation  within  r,  and  r  =  lx|.  Up  is  an  irrotational  velocity 

2 

field  added  to  satisfy  the  potential  boundary  conditiorf.  Up  =  V<J>  where  V<f>  =0 
and  u.n  =  0  on  solid  boundaries  while  u.n  =  U  at  the  inlet,  n  is  the  normal 
unit  vector.  For  the  confined  shear  layer,  the  boundary  condition  at  x  =  0  is: 
u  =  U1  for  y  is  greater  than  0  and  u  =  U2  for  y  less  than  0,  while  at  y  =  0 
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there  is  a  vortex  sheet  of  strength  AU  =  U1-U2. 

In  this  case,  we  use  Rankine  vortex  elements,  i.e.,  the  vorticity  of  an 
element  is  constant  within  the  core  and  zero  outside,  f(r)  =  1/jt  for  r  <  6  and 
f(r)  =  0  fo»"  ^  >  6.  Correspondingly,  k(r)  =  r^/2rt  for  r  <  6  and  k  =  1  for  r  > 
6.  Moreover,  the  potential  velocity  field  is  obtained  by  conformal 
transformation.  Thus,  the  physical  plane  is  mapped  onto  the  upper  half  plane 
and  image  vortices  are  used  to  satisfy  the  potential  boundary  conditions.  The 
form  of  the  mapping  function  for  the  confined  shear  layer  is  given  by  Ghoniem 
and  Ng  (1986). 

The  motion  of  the  vortex  elements  must  be  constructed  such  that  the 
vorticity  field  satisfied  Equation  (1).  This  is  accomplished  by  solving  this 
equation  in  two  fractional  steps: 

convection:  -fr  +  u-  7(u  =  0  (7) 

o  t 

diffusion:  77  =  ~~72(»  (8) 

d  t  Re 


In  the  first  step,  the  convective  transport  of  vorticity  is  implemented  in 
terms  of  the  Lagrangian  displacement  of  the  vortex  elements  using  the  current 
velocity  field  computed  from  Equation  (6).  In  the  second  step,  the  solution  of 
the  diffusion  equation  is  simulated  stochastically  by  the  random  walk 
displacement  of  the  vortex  elements  according  to  the  appropriate  population. 
Thus, 

Xi  (t  +  At)  =  Xi(t)  +  ^  s(Xik)  At  +  ui  (9) 

for  i  =  1.2, ...,h,  where  ^  is  a  k-th  order  time- integration  scheme  and  iu  is  a 
two-dimensional  Gaussian  random  variable  with  zero  mean  and  standard  variation 
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/2  a  t  /  Re\  For  more  details,  see  Ghoniem  and  Ng  (1986),  Ghoniem  and  Gagnon 
(1986) . 

The  no-slip  boundary  condition  at  the  walls  is  satisfied  by  generating  new 

vortex  elements  to  cancel  the  induced  velocity  by  the  vorticity  field.  Here, 

we  generate  vorticity  only  at  the  point  of  separation,  i.e.,  at  the  tip  of  the 

splitter  plate,  since  the  growth  of  the  boundary  layers  along  the  channel  walls 

at  these  high  Reynolds  numbers  is  small.  At  each  time  step,  the  new  vorticity 

AT  =  -  AUU  At, where  Um  =  (U1  +  U2)/2,  is  consigned  to  No  elements  of  strength 
m 

AT/No  and  added  to  the  field  at  points  Ax  =  UmAt/No  apart  downstream  of  x  =  0 
(Figure  1). 

The  effect  of  the  numerical  parameters  on  the  accuracy  of  the  results  was 
investigated  by  Ghoniem  and  Ng  (1987).  Their  results  emphasized  the  importance 
of  using  a  high  order  time- integration  scheme  with  k=2  to  avoid  excessive 
numerical  diffusion  in  the  vorticity  field.  The  value  of  No  =  6  was  also  found 
to  be  appropriate,  in  order  to  obtain  well-defined  eddy  structures  after  the 
rollup  and  the  first  two  pairings.  The  second  pairing  is  accomplished  within 
the  domain  of  0  <  x  <  6,  therefore,  the  computational  domain  was  limited  to 
Xmax  =  6.  Downstream  of  Xmax,  the  vorticity  was  deleted.  Varying  Xmax  showed 
that  the  effect  of  deleting  the  vortex  elements  propagates  about  one  channel 
height  upstream,  hence,  the  results  are  accurate  only  for  x  <.  5. 

Flow  Structure 

To  address  the  capabilities  of  the  vortex  method,  from  flow  visualization 
point  of  view,  the  results  of  some  sample  simulations  are  presented  in  this 
section. 

Results  of  a  typical  simulation,  presented  in  terms  of  the  velocity  and 
location  of  all  vortex  elements  used  in  the  computations,  are  shown  in  Figures 
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2,  3  and  4  for  the  cases  of  Re  =  24000,  Re  =  4000,  and  Re  =  1000,  respectively. 
Each  vortex  element  is  depicted  by  a  point,  while  its  velocity  relative  to  the 
mean  velocity  is  represented  by  a  line  vector  starting  at  the  center  of  the 
vortex  element.  The  velocity  ratio  across  the  layer  at  the  inlet  is  U2/U1  = 
1/3. 

Results  show  the  formation  of  large  vortex  eddies  by  the  rollup  of  the 
vorticity  layer  that  emanates  at  the  splitter  plate  and  the  subsequent  pairings 
of  these  eddies  into  larger  structures.  The  rollup  of  the  shear  layer  was 
investigated  in  Ghoniem  and  Ng  (1986)  by  analyzing  results  at  a  wide  range  of 
the  Reynolds  number  and  at  different  boundary  conditions.  Their  analysis  show 
that:  (1)  the  rollup  is  due  to  the  growth  of  perturbations  by  the  Kelvin- 
Helmholtz  instability  mechanism,  and  the  shedding  frequency  corresponds  to  the 
most  unstable  frequency  predicted  from  the  linear  stability  analysis  of  a 
spatially  growing  layer*,  (2)  pairing,  which  is  associated  with  the  local  sub¬ 
harmonic  perturbations,  results  in  a  step-wise  increase  in  the  size  of  the 
vorticity  layer  as  two  eddies  merge',  (3)  the  two  sources  of  the  subharmonic 
perturbations  are  the  downward  motion  of  the  layer  and  the  monotonic  growth  in 
the  size  of  the  eddies  downstream',  (4)  the  intrinsic  dynamics  of  the  instabil¬ 
ity  is  not  strongly  affected  by  the  value  of  the  Reynolds  number,  except  that 
at  the  low  Reynolds  number  the  eddies  are  slightly  larger  due  to  the  dispersion 
of  vorticity  by  diffusion,  and  (5)  the  computed  velocity  statistics  show  good 
agreements  with  experimental  data,  indicating  that  the  fundamental  mechanisms 
of  the  shear  layer  are  two-dimensional  and,  hence,  the  numerical  scheme  is 
capable  of  predicting  the  large  scale  features  accurately. 

3.  IMPROVED  VORIEX  MEIHODS 

Despite  the  encouraging  results  obtained  by  vortex  methods  in  exhibiting 
the  largt  scale  features  of  the  unstable  shear  layer,  it  has  been  recently 
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found  that  using  a  fixed  number  of  vortex  elements  can  lead  to  large  errors  in 
such  approaches  (Ghoniem  et  al.  1987a).  These  errors  are  mainly  caused  due  to 
the  strong  strain  field  that  develops  and  acts  to  distort  the  original  vortic- 
ity  contours.  Analysis  of  the  convergence  properties  of  inviscid  vortex 
methods  has  revealed  that  the  ratio  of  the  core  radius  to  the  separation 
between  vortex  elements  is  a  major  factor  in  governing  the  accuracy  of  vortex 
computations  (Chorin,  1973*.  Hald,  1979’,  and  Beale  and  Majda,  1982).  Another 
factor  influencing  the  accuracy  are  the  form  of  the  core  function  and  the 
scheme  of  discretizing  the  initial  vorticity. 

Figures  2,  3,  and  4  clearly  show  that  as  the  magnitude  of  the  strain  rate 
increases,  due  to  the  rollup  of  the  unsteady  shear  layer,  the  separation 
distance  within  the  vortices  increases.  Considering  the  fact  that  these 
calculations  are  based  on  a  fixed  number  of  elements  with  invariant  core  sizes, 
this  may  result  in  large  errors.  Large  local  strains  are  created,  which  would 
lead  to  strong  plane  stretch  of  the  vorticity  and  would  cause  the  initial 
circular  patch  of  vorticitiy  to  deform  into  elliptical,  or  more  complex,  shapes 
with  their  major  axes  aligned  with  the  principle  direction  of  the  strain  field. 
If  the  number  of  computational  elements  that  maintain  their  original  shapes  is 
used,  they  would  not  be  able  to  accommodate  these  severe  changes,  and  this 
problem  may  threaten  the  long-time  behavior  of  the  method.  Using  a  large 
number  of  elements  to  discretize  the  initial  patch  of  vorticity,  these  elements 
will  naturally  migrate  in  the  proper  direction,  Le.,  direction  of  high  strain, 
and  complex  formation  will  arise  without  the  need  to  modify  the  cores  of 
individual  elements.  However,  that  may  necessitate  the  use  of  a  prohibitively 
large  number  of  elements  initially,  or  doing  a  series  of  runs  while  increasing 
the  number  of  elements  to  guarantee  the  presence  of  enough  elements  after  the 
development  of  the  large  strain  field  until  a  convergent  solution  is  obtained. 
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Id  another  work,  Ghoniem  et  al.  (1987b)  show  that  using  extra  elements  in 
the  direction  of  local  stretch  is  a  good  way  of  preserving  the  accuracy.  When 
the  distance  between  two  neighboring  elements  in  the  direction  of  principle 
strain  exceeds  a  minimum  value,  a  new  element  is  placed  between  the  two.  This 
procedure  should  be  done  in  such  a  way  that  the  total  circulation  is  conserved 
and  the  vorticity  distribution  is  not  strongly  perturbed.  In  the  simplest 
form,  one  may  replace  the  two  vortex  elements  on  a  line  by  three  elements  along 
the  same  line  (Figure  5a).  This,  obviously,  might  result  in  an  excessive 
number  of  computational  elements  within  the  domain  of  interest.  To  remedy 
this,  one  may  similarly  extract  elements  and  combine  the  two  neighboring 
elements  if  their  separating  distance  decreases  below  a  critical  minimum  value 
(Figure  5b).  Ike  critical  values  for  the  distance  between  the  neighboring 
elements  to  require  the  injection  of  a  new  element,  or  to  combine  the  two 
elements  into  a  single  one,  determines  the  accuracy  of  the  calculations  and  is 
obviously  dependent  on  the  computational  capabilities.  Moreover,  this 
algorithm  requires  more  elaborate  programming  than  conventional  vortex  methods, 
since  it  uses  information  about  neighboring  elements  in  the  principle  direction 
of  strain  (Ghoniem  et  al.,  1987a). 

To  implement  the  scheme,  we  revisited  the  spa  ti  al  ly-  evol  ving  flow  and 
distributed  the  vorticity  at  the  inlet  among  five  layers  of  vortex  elements. 
This  distribution  was  represented  by  a  Gaussian  curve  with  a  spread  2a: 

G(x)  =  exp  (-  ^5-)  (10) 

/no  a 

which  corresponds  to  an  Error  function  distribution  for  the  velocity. 

To  ensure  the  accuracy  of  discretization  of  the  vorticity  among  the  five 
layers,  it  was  shown  that  either  a  pointwise  discretization,  uk  =  Q(X^),  or  an 
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ares  average  approach,  uk  =  J  Q(x)  dx  (where  Xi  are  the  centers  of  a  square 
mesh  of  size  h)  to  distribute  the  vorticity  of  the  shear  layer  among  vortex 
elements  produce#  1  arge  error.  Instead,  the  following  scheme  was  used: 

N 

Q(Xi)  =  u.h2  f6(Xi  -  Xj)  (11) 

j=l 


where  Ti  =  ckI*  is  the  circulation  of  an  element  of  strength  wi,  6  is  the 
core  radius,  and  fc  is  the  core  function.  The  function  f«>(r)  i*  written  as 


f&(x)  =  K  f(r/6) 
0  6Z 


(12) 


where 


r2  =  x2  +  y2 


(13) 


The  function  fg  is  radially  symmetric,  i.e.,  J  fgdx=l  and  is  a  fast  decaying 
function  of  tJ 6  so  that  most  of  the  vorticity  is  concentrated  within  r  <  6.  In 
the  calculations  presented  below,  we  used  a  second  order  Gaussian  core,  as 
suggested  by  Beale  and  Majda  (1985). 

i  2 

f ( r)  =  -  e~r  (14) 

TT 

In  Figure  6,  the  location  and  velocity  of  all  the  vortex  elements  are 
presented.  In  this  figure,  the  time  sequence  of  the  vorticity  plot  exhibits  the 
large  scale  features  of  the  flow  clearly.  In  this  calculation,  the  effects  of 
Reynolds  number  is  not  taken  into  account  (i.e.,  inviscid  calculations)  and  the 
values  of  a  and  P  (Figure  5)  were  chosen  to  be  low  near  the  inlet,  and  to 
increase  linearly  with  the  streamwise  direction  x.  This  is  described  by: 

a  =  ao  +  Aa  x  (15) 

p  =  po  +  A0  x  (16) 

For  the  calculations  in  Figure  6,  aQ  =2,  po  =  1,  Aa  =  3  and  AP  =  2. 
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There  are  no  explicit  perturbations  added  onto  this  flow,  and  similar  to  the 
cases  of  Figures  2  through  4,  the  numerical  truncation  errors  are  the  main 
reason  for  perturbing  the  flow  and  causing  the  rollup  of  the  unsteady  shear 
layer.  There  are  approximately  600  vortex  elements  within  the  domain  and,  as 
can  be  seen  in  Figure  6,  the  initial  stages  of  the  development  of  the  unsteady 
shear  layer  is  not  clearly  resolved  by  the  selected  values  of  a  and  p.  A 
finer  resolution  calculation  was  performed  by  decreasing  the  values  of  Aa  and 
Ap  to  2.5  and  1.5,  respectively,  while  keeping  the  values  of  a0  and  P0 
unchanged.  The  results  of  this  calculation  are  presented  in  Figure  7,  and  a 
comparison  between  this  figure  and  Figure  6  indicates  that  the  structure  of 
large  scales  of  the  flow  is  much  more  clearly  recognizable.  In  this  case, 
there  are  about  1000  vortex  elements  within  the  computational  domain  and  the 
simulations  required  about  45  minutes  of  cpu  time  on  the  Wright-Patterson  CRAY- 
XMP  to  reach  a  normalized  computational  time  of  t  =  36.  Higher  resolution 
simulation  with  aQ  =  1.5,  po  =  0.5,  Aa  =  2.5  and  Ap  =  1.5  shows  the  large  scale 
feature  of  the  flow  even  more  clearly.  This  is  presented  in  Figure  8,  where 
the  location  and  the  relative  velocities  of  approximately  2000  elements  are 
shown.  Future  calculations  with  better  resolution  (smaller  values  of  a  and  p) 
are  recommended  to  test  the  dependency  of  the  final  results  on  the  process  of 
vortex  injections  and  recombinations.  Harmonic  perturbations  would  speed  up 
the  processes  of  the  vortex  rollup  and  pairing,  as  shown  in  Figure  9.  In  this 
figure,  the  response  of  the  layer,  with  the  same  hydrodynamical  parameters  as 
that  in  Figure  6,  to  the  application  of  the  most  unstable  frequency  with  low 
amplitude  forcing  on  the  upstream  side  of  the  shear  layer  is  presented.  The 
vorticity  plots  of  Figures  6  and  7  indicate  that  the  effect  of  harmonic  forcing 
is  to  increase  the  mixing  rate  in  the  initial  stages  of  mixing  immediately 
after  the  splitter  plate,  where  the  resonant  eddy  is  forming.  This  region  is 
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marked  by  a  row  of  organized  structures,  which  are  created  approximately  at 
equal  wavelengths  from  each  other  due  to  the  application  of  harmonic  forcing 
(Figure  9)  (Brow  and  and  Ho,  19*  .*,  Oster  and  Wygnanski,  1982).  At  this  region, 
the  increased  mixing  rate  results  in  the  augmentation  of  the  chemical  reactions 
(if  two  reactants  were  introduced  through  the  two  streams  of  the  layer).  The 
results  of  these  calculations  are  consistent  with  the  previous  results  of  Givi 
and  J ou  (1988),  who  showed  that  the  harmonic  forcing  plays  an  important  role  in 
the  active  control  of  the  combustion  in  the  mixing  region  of  the  layer  if  the 
reaction  is  assumed  to  be  under  equilibrium.  What  happens  for  nonequilibrium 
conditions,  when  the  increased  mixing  might  result  in  decreased  reactions  (Givi 
et  al.,  1987),  cannot  be  investigated  in  the  present  calculations  and  is  the 
subject  of  our  future  inve stiga tion. 

As  mentioned  before,  a  particular  advantage  of  the  improved  vortex  method 
is  in  its  capability  in  resolving  the  regions  of  high  scalar  gradients.  This 
superiority  can  be  recognized  by  a  comparison  between  the  vorticity  plots  in 
Figures  7  through  9  and  those  in  Figures  2  through  4.  This  comparison  show  s 
that  the  application  of  the  previously  used  vortex  methods  cannot  resolve  the 
of  high  strain,  as  the  vortices  tend  to  evacuate  these  regions  (i.e.,  braids  of 
the  vortices  in  Figures  2  through  4).  The  improved  vortex  method,  however,  has 
the  advantage  that  the  injection  procedure  (Figure  5a)  prevents  the  vortices 
from  moving  away  from  the  braids  and,  therefore,  allows  for  the  accurate 
simulations  of  those  regions  of  the  flow  field  under  large  local  stretch.  This 
feature  of  the  improved  vortex  method  would  be  particularly  useful  in  diffusion 
flame  calculations  where  the  chemical  nonequilibrium  effects  are  most  important 
at  regions  of  high  dissipation  (strain)  rates  (Givi  et  al.,  1987*,  Peters,  1984*, 
Tsuj  i,  1982,  andLiew  et  al.,  1984). 
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The  Transport  Element  Method 

The  same  philosophy  adopted  in  the  vortex  methods  for  the  hydrodynamic 
simulation  can  be  used  for  the  transport  of  scalar  quantities.  In  such 
approach,  the  regions  of  high  scalar  gradients  may  be  discretized  by  transport 
elements  (particles)  carrying  a  certain  amount  of  the  gradient  of  the  scalar. 
This  is  equivalent  to  replacing  u  and  Re  in  Equation  (1)  by  Q  =  VT  and  Pe, 
where  1  is  a  generalized  scalar  (temperature  or  concentration)  and  Pe  is  the 
Peclet  number.  However,  as  will  be  seen  in  this  section,  contrary  to 
vorticity,  the  scalar  gradients  are  not  conserved  along  particle  path,  and 
their  strengths  should  be  modified  according  to  the  local  stretch  and  tilting 
of  the  material  element.  To  explain  this  approach,  the  transport  of  the 
conserved  scalar  quantity,  T,  is  considered 

a-  vT  -  v2t  (17) 

using  the  fractional  stepping  method,  the  same  as  that  in  Equations  (7)  and 


(8),  the  solution  of  Equation  (17)  is  equivalent  to: 


Convection  TT  +  -B'  VT  =  0 
d  t 


Diffusion  |f  =  V*T 


The  transport  element  method  is  implement  by  taking  the  gradient  of  the 
Equations  (18)  and  (19) 


dO 

—  +  jg.yQ+Q.Vu+  Qx(Vxfi)  +  ux 


(VxQ) 


0 

0 


The  last  term  in  Equation  (20)  is  zero.  Since  Q  is  the  gradient  of  a  scalar 
quantity.  For  a  two-dimensional  flow,  Q  can  be  written  as  Q  =  (p,  q),  where  p 
and  q  are 

«•»  -  <g  ■  g> 

In  terms  of  q  and  p.  Equations  (20)  and  (21)  become: 


£<l 

at 


n-Vq 


(23) 


fljE 


+  u'Vp 


3u  dv 
q  3y  ^  3y 


(24) 


19  _  J_  ^2 
at  Pe  v 


q 


(25) 


»  -  k  *  ,26) 

The  RHS  of  Equations  (23)  and  (24)  indicate  that  the  strength  of  the  scalar 
gradients  (i.e.,  the  magnitude  of  p  and  q  within  each  element)  would  change 
along  the  particle  path.  Therefore,  the  convective  step  (Equations  23  and  24) 
will  be  completed  in  two  fractional  steps: 


Pure 

convecti on 


+ 

at 


£•  Vq  =  0 


(27) 


Source 
adj  ustment 


3  q  du  dv 

dt  **  3y  ^  3x 


(28) 


A  similar  procedure  would  also  follow  for  p. 

In  the  first  fractional  step  (Equation  27),  the  location  of  the  gradient 
element  changes  in  aLagrangian  sense,  and  in  the  second  fractional  step  the 
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strength  of  the  gradient  is  changed  to  satisfy  Equation  (28).  The  magnitude  of 
the  gradient  at  a  location  x  is,  then,  computed  from: 

N 

a(x,t)  =  ^  Qi(t)h2f6  (x  -  xA)  (29) 

where  f  is  the  core  function  for  the  transport  elements  and  can  be  chosen  to  be 
different  from  the  vorticity  elements.  In  this  work,  we  used  the  same  core 
functions  (second  order  Gaussian). 

Once  the  location  and  the  magnitude  of  the  strength  of  the  scalar 
gradients  (p,q)  are  known,  the  magnitude  of  the  scalar  (T)  can  be  easily 
obtained  by  using  the  Green's  function  approach  as  developed  by  Anderson 
(1983).  According  to  this  method,  the  solution  of  the  equation 

V21  =  V  -  Q  (30) 

can  be  written  as 

T  =  J  (7-fi)  Gdx  (31) 

where 

G  =  -  ~  In  r  (32) 

is  the  Green  function  of  the  Poisson  equation.  Integration  Equation  (31)  by 
parts  results  in: 

T  =  J  (fi*VG)  dx  (33) 

Using  Eq.  (29)  for  fi,  we  have: 

N 

Kx.t)  =  )  fii(t)*7Gc  (£  -  jSi(t))h2  < 34) 

i=l 
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where 


G6  ‘  ^  k<l>  (35) 

r 

and  6  is  the  core  function. 

This  expression  is  convenient  when  used  in  connection  with  the  vortex 
method  because  the  same  elements  used  in  the  transport  of  the  vorticity  may  be 
used  for  the  scalar  gradient  transport  (when  Re  =  Pe  and  the  same  initial 
boundary  conditions).  Also,  the  same  injection-recombination  procedure  as  used 
in  the  improved-vortex  calculations,  may  be  used  in  the  gradient  transport 
simulations  to  accurately  simulate  the  regions  of  high  scalar  gradients.  The 
need  for  the  insertion  is  more  apparent  here  since  the  magnitude  of  the 
gradient  increases  where  the  strain  field  is  high. 

To  implement  the  procedure,  the  transport  of  a  conserved  scalar  variable 
with  a  normalizedvalue  of  T  =  0  in  the  high  speed  stream  and  T=  1  in  the  low 
speed  stream  is  considered.  The  magnitude  of  the  diffusivity  of  this  scalar 
quantity  is  assumed  zero  so  that  the  same  vortex  elements  can  be  used  for  the 
transport  of  the  gradient  Q.  Calculations  were  performed  for  the  same 
hydro  dynamical  conditions  as  that  presented  in  Figure  6.  Therefore,  the  vortex 
plots  in  this  figure  would  also  exhibit  the  gradients  plots  in  that  the 
locations  of  the  transport  elements  and  those  of  the  vortices  coincide.  Of 
course,  the  strengths  of  the  computational  elements  would  be  different  in  the 
two  cases. 

The  instantaneous  values  of  the  scalar  variable  T,  calculated  from 
Equation  (34)  are  statistically  analyzed  for  the  purpose  of  comparison  with 
experimental  data.  This  analysis  was  performed  by  averaging  the  instantaneous 
values  at  every  time  step  during  the  period  tl  <  t  <  t2.  The  vlaue  of  tl  was 
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selected  as  twice  the  value  of  the  residence  time  tl  =  2xmax/0m,  so  that  the 
effects  of  initial  transients  are  washed  out  of  the  domain.  The  period  of  AT  = 
tl-t2  was  selected  as  four  times  the  value  of  the  residence  time  (i.e.,  AT  = 
4Xmax/Um),  so  that  it  is  lsrge  enough  for  the  purpose  of  statistical  analysis. 
During  this  period  an  ensemble  of  160  data  points  was  gathered  that  resulted 
In  almost  identical  statistics  to  those  with  an  ensemble  of  80  sample  data 
points.  Therefore,  higher  ensemble  sizes  were  not  considered  to  minimize  the 
computational  costs. 

The  mean  and  rms  values  of  the  scalar  variable  T  are  presented  in  Figures 
10  and  11  versus  the  reduced  coordinate,  =  (Y-Yo) / (X-Xo) ,  where  Yo  is 
measured  at  T  =  0.5  (-  indicate*  the  ensemble  average)  and  Xo  is  the  virtual 

origin  of  the  mixing  layer  based  on  the  mean  scalar  profile  (in  this 
simulation,  Xo  =  0).  In  these  figures,  the  solid  lines  are  the  results  of 
simulations  and  the  data  points  are  obtained  from  the  experimental  measurements 
of  Masutani  and  Bowman  (1986)  for  a  dilute  nonreacting  mixing  layer  with  the 
same  velocity  ratio.  It  is  evident  from  these  two  figures  that  both  the  mean 
and  the  rms  of  the  passive  scalar  are  accurately  predicted  across  the  width  of 
the  shear  layer.  It  is  particularly  encouraging  to  note  that  both  these 
profiles  are  in  better  agreement  with  data  than  those  obtained  previously  by 
Givi  et  al.  (1985)  utilizing  a  K-e  turbulent  model  and  a  gradient  diffusion 
model  of  turbulent  transpor  t  of  the  scalar  mean,  moments  and  PDF  (probability 
density  function).  In  particular,  the  K-e  model  was  not  capable  of  accurately 
predicting  the  behavior  of  the  scalar  field  near  the  high  speed  stream  where 
the  effects  of  large  scale  structures  are  important.  The  effects  of 
interm i t tency,  which  are  not  included  in  the  K— E  calculations,  are  explicitly 
accounted  for  in  the  present  simulations  (Figures  6  through  9)  and  allow  for 
better  comparison  with  data  near  the  high  speed  stream.  Moreover,  the 
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simulations  indicate  the  presence  of  two  local  maxima  in  the  rms  profiles  that 
correspond  approximately  to  the  locations  where  the  gradient  of  the  mean  valoe 
is  highest,  as  also  indicated  in  the  results  of  Masutani  and  Bowman  (1986)  and 
Batt  (1977).  In  view  of  these  comparisons  with  data,  it  is  clear,  in 
accordance  with  the  findings  of  Broadwell  and  Briedenthal  (1982),  that  the 
intermittency  caused  by  the  large  coherent  structures  contribute  greatly  to  the 
statistics  of  this  unsteady  flow. 

The  probability  density  functions  (PDF's)  of  the  scalar  variable  was  also 
constructed  from  the  160  ensemble  data  points.  Ihe  cross  stream  variation  of 
the  PDF's  are  presented  in  Figure  12  for  the  same  streamwise  location  as  those 
in  Figures  10  and  11.  In  this  figure,  T  represents  the  instantaneous  scalar 
values  and  the  PDF  is  defined  such  that 

P(T)  dT  =  Prob  {x  <  T  <  T  +  dT)  (36) 

In  this  figure,  the  delta  function  at  T  =  0  indicates  the  fluid  from  the  high 
speed  stream  and  the  delta  function  at  T  =  1  represents  the  low  speed  fluid  con¬ 
centration.  It  is  also  observed  that,  as  the  layer  is  traversed,  the  height  of 
the  PDF's  within  the  layer  represents  the  fluid  in  the  mixing  core.  The  results 
presented  in  this  figure  are  similar  to  those  observed  experimentally  in  the 
pretransitional  regions  of  the  layer  (Koochesf ahani,  1984).  However,  they  do 
not  agree  with  the  recent  studies  of  Givi  and  McMurtry  (1987)  and  Givi  and  Jou 
(1988),  who  showed  that  the  mechanisms  of  the  vortex  rollup,  and  the  pairing 
process  in  particular,  contribute  greatly  to  the  generation  of  a  third  peak  (a 
third  delta  function)  at  a  mixture  concentration  whose  value  depends  on  the 
velocity  difference  across  the  layer.  There  are  two  main  reasons  for  this  dis¬ 
agreement.  First,  in  the  present  calculations,  only  the  convective  transport 
of  the  scalar  is  considered  and  the  contributions  of  the  diffusive  transport 
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are  taken  into  account  (i.e.,  infinite  Peclet  number).  The  effects  of 
diffusivity  are  to  expand  the  core  of  the  transport  element  and.  as  a  result, 
to  mix  the  fluid  better  within  the  mixing  region.  That  may  contribute  to  the 
generation  of  the  predominant  peak  of  the  PDF's.  Second,  the  number  of  vortex 
(and  transport)  elements  employed  in  the  present  calculations  may  not  be  large 
enough  to  accurately  capture  the  PDF's,  which  are  "sensitive"  statistically 
even  though  less  sensitive  quantities,  such  as  the  mean  and  the  rms,  are 
accurately  predicted. 

Future  calculations  with  larger  number  of  computational  elements  will  be  a 
logical  continuation  of  the  present  work. 

4 .  CONCLUSION 

A  Lagrangian  numerical  scheme  based  on  an  improved  v or tex- tr anspor t 
element  method  has  been  employed  to  investigate  the  mechanism  of  mixing  in  a 
spatially  developing,  two-dimensional  shear  layer  configuration.  In  this 
scheme,  Lagrangian  elements  are  used  to  transport  vorticity  and  also  the 
gradients  of  the  scalar  variables.  A  particular  advantage  of  this  method  is 
that  it  optimizes  the  computational  efforts  by  distributing  the  vortex  and  the 
transport  elements  around  the  regions  of  high  vorticity  and  high  gradients. 

An  injection  procedure  is  employed  to  prevent  the  computational  elements 
from  evacuating  the  regions  of  high  local  stretch.  The  implementation  of  this 
procedure  makes  this  improved  method  more  attractive  than  the  conventional 
vortex  methods  in  that  regions  of  large  dissipation  rates  (i.e.,  the  braids  of 
the  vortices)  can  be  more  accurately  simulated.  Ihe  coherent  large  structures 
of  the  flow  can  be  nicely  exhibited  by  "fine  tuning"  of  this  injection 
procedure. 

The  resulting  scheme  was  utilized  to  simulate  the  convective  transport  of 
a  passive  scalar  quantity  within  the  mixing  region  of  the  layer.  The 
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statistical  analysis  of  the  instantaneous  results  showed  favorable  agreements 
with  experimental  data  for  the  first  two  moments  (mean  and  ms)  of  the  scalar. 
However,  it  seems  that  to  make  a  better  comparison  with  measured  data  for  the 
PDF's  of  the  scalar  variable,  simulations  should  be  performed  with  larger 
number  of  transport  elements. 

Harmonic  perturbation  of  the  velocities  at  the  inlet  of  the  layer  was 
shown  to  have  a  significant  effect  in  accelerating  the  process  of  the  vortex 
rollup  and  creating  large  organized  structures  within  the  mixing  regions  of  the 
flow.  This  resulted  in  the  convolution  of  the  interphase  surfaces  and  the 
improvement  of  mixing  in  the  region  immediately  behind  the  splitter  partition 
plate.  This  would  be  very  important  in  controlling  the  outcome  of  the  chemical 
reactions  in  combusting  shear  layer  configuration,  which  would  be  a  logical 
extention  of  this  work. 

Calculations  with  higher  number  of  transport  elements  are  recommended  also 
for  better  flow  visualizations  and  to  check  the  sensitivity  of  the  statistical 
variables  on  the  total  number  of  Lagrangian  elements. 
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Figure  2.  Vorticity  field  at  Re  -  24000  and  U2/Ui  -  1/3  by  vortex  method.  The  splitter 
plate  is  located  on  the  left  side  of  the  figure. 
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Figure  3.  Vorticity  field  at  Re  ”  4000  and  02^1  “  1/3  by  vortex  method.  The  splitter  plate 
is  located  on  the  left  aide  of  the  figure. 
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Figure  A.  Vorticity  field  at  Re  ”  1000  and  02^1  *"  1/3  by  vortex  method.  The  splitter  plate 
la  located  on  the  left  side  of  the  figure. 
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Figure  5.  Injection  tod  recombination  mechanists  in  improved  vortex  method. 
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Figure  6.  Vorticity  field  by  improved  vortex  method.  Inviacid  calculations.  The  splitter  plat 
is  located  on  the  left  side  of  the  figure.  an  m  2,  ■  1,  a<*  ■  3,  &6  ■  2. 
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Vorticlty  field  by  improved  vortex  method.  Inviscid  calculations.  The  splitter  plate 
is  located  on  the  left  side  of  the  figure.  a0  ■  1.5,  B0  -  0.5,  Aa  •  2.5,  AB  ■  1.5. 
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Figure  9.  Vorticity  field  by  improved  vortex  method.  Inviscid  calculations.  Harmonically  forced 
layer  at  the  same  conditions  as  in  Figure  6.  The  splitter  plate  is  located  on  the  left 
side  of  the  figure. 
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Figure  11.  Ras  scalar  profile  vs.  the  croBs-stream 
coordinate.  ,  present  simulations; 

0  ,  A  ,  □  ,  experimental  data  obtained 
by  Masutani  and  Bowman  (1986)  at  different 
stream-nrlse  locations. 
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Figure  12.  PDF's  scslsr  vs.  the  scslsr  space  and  the  cross-stress  coordinate. 


